Set Up Steps

library(censusapi)
library(tidycensus)
library(tidyverse)
library(sf)
library(geojsonsf)
library(mapview)
library(dplyr)
library(plotly)
library(tigris)
library(readxl)
library(leaflet)
library(RColorBrewer)
library(sp)
library(usmap)

mapviewOptions(
  basemaps = "CartoDB.Positron",
  #vector.palette = colorRampPalette(cols)
)

options(
  tigris_class = "sf",
  tigris_use_cache = TRUE
)

Sys.setenv(CENSUS_KEY="0c313bd613a7281ae62c2fbb004d156d647e9c94")
projection <- "+proj=utm +zone=10 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=ft +no_defs"
# Obtain the census variable lookup data
# acs_vars <-
#   listCensusMetadata(
#     name = "2018/acs/acs5",
#     type = "variables"
#   )

# Save an .rds version
# saveRDS(acs_vars, file = "/Users/armellecoutant/Documents/GitHub/218z/Armelle_Coutant/sj_outreach/censusData2018_acs_acs5.rds")

# Obtain the saved census data
acs_vars = readRDS("/Users/armellecoutant/Documents/GitHub/218z/Armelle_Coutant/sj_outreach/censusData2018_acs_acs5.rds")
# Get Santa Clara block groups
sc_blockgroups <-
  block_groups("CA","Santa Clara", cb=F, progress_bar=F) %>% 
  st_transform('+proj=longlat +datum=WGS84')

# Get San Jose geometry
sj_geom <- places("CA", cb = F, progress_bar = FALSE) %>% 
  filter(NAME == "San Jose") %>% 
  st_transform('+proj=longlat +datum=WGS84')

# Filter to Santa Clara block groups within San Jose geometry
sj_blockgroups <- sc_blockgroups[st_contains(sj_geom, sc_blockgroups %>% st_centroid())[[1]],]

Vulnerability Mapping

These vulnerability maps are intended to determine the San Jose block groups with limited access to information and resources. The following factors were selected for the layered map: - Households without Internet Access - Households without Computers - Households Below 80% AMI - K-12 Enrollment - Lack of Small Cell Coverage

Households Without Internet Access

sc_internet <- 
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B28002)"
  ) %>%
  mutate(
    blockgroup = paste0(state,county,tract,block_group)
  ) %>%
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
  dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(
    key = "variable", 
    value = "estimate", 
    -blockgroup
  ) %>% 
  mutate(
    label = acs_vars$label[match(variable,acs_vars$name)]
  ) %>% 
  dplyr::select(-variable) %>%
  separate(label, into = c(NA, NA, "subscription", "type", "additional"), sep = "!!") %>% 
  filter(is.na(subscription) | subscription == "No Internet access") %>%
  mutate(
    subscription = replace_na(subscription, "total_num")
  ) %>%
  dplyr::select(blockgroup, estimate, subscription) %>%
  spread(key = subscription, value = estimate) %>%
  mutate(
    percent_no_internet = (`No Internet access`*100 / total_num) %>% round(1)
  )
sj_blockgroups_internet <- 
  sc_internet %>% 
  filter(blockgroup %in% sj_blockgroups$GEOID) %>% 
  left_join(sj_blockgroups, by = c("blockgroup" = "GEOID")) %>% 
  st_as_sf() %>%
  st_set_crs(4326) %>%
  st_transform(projection)

sj_blockgroups_internet[is.na(sj_blockgroups_internet)] = 0

total_no_internet_sj = sum(sj_blockgroups_internet$`No Internet access`)
total_sj <- sum(sj_blockgroups_internet$total_num)
perc_no_internet_sj <- total_no_internet_sj*100/total_sj

This is the total number of households in San Jose without internet:

print(total_no_internet_sj)
## [1] 24834

This is the percent of households in San Jose without internet:

print(perc_no_internet_sj)
## [1] 8.052816

This map highlights San Jose block groups by count of households without internet.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_blockgroups_internet %>% 
      pull(`No Internet access`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_blockgroups_internet %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`No Internet access`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    label = paste0(sj_blockgroups_internet$`No Internet access`)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_blockgroups_internet$`No Internet access`,
    title = 'Households<br>Without Internet'
  )

This is a comparable map, except that households without internet are represented as a percent of the total number of households within the block group.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_blockgroups_internet %>% 
      pull(`percent_no_internet`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_blockgroups_internet %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`percent_no_internet`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    label = paste0(sj_blockgroups_internet$`percent_no_internet`)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_blockgroups_internet$`percent_no_internet`,
    title = 'Percent<br>Without Internet'
  )

Households Without a Computer

sc_computer <- getCensus(
  name = "acs/acs5",
  vintage = 2018,
  region = "block group:*", 
  regionin = "state:06+county:085",
  vars = "group(B28003)"
  ) %>%
  mutate(blockgroup = paste0(state,county,tract,block_group)) %>%
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>%
  dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>% 
  gather(key = "variable", value = "estimate", -blockgroup) %>% 
  mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>% 
  dplyr::select(-variable) %>%
  separate(label, into = c(NA, NA, "computer", "internet"), sep = "!!") %>% 
  filter(is.na(computer) | computer == "No computer") %>%
  mutate(computer = replace_na(computer, "total_num")) %>%
  dplyr::select(blockgroup, estimate, computer) %>%
  spread(key = computer, value = estimate) %>%
  mutate(percent_no_computer = (`No computer`*100 / total_num) %>% round(1))
sj_blockgroups_computer <- 
  sc_computer %>% 
  filter(blockgroup %in% sj_blockgroups$GEOID) %>% 
  left_join(sj_blockgroups, by = c("blockgroup" = "GEOID")) %>% 
  st_as_sf() %>%
  st_set_crs(4326) %>%
  st_transform(projection)

sj_blockgroups_computer[is.na(sj_blockgroups_computer)] = 0

total_no_computer_sj = sum(sj_blockgroups_computer$`No computer`)
total_sj_comp <- sum(sj_blockgroups_computer$total_num)
perc_no_computer_sj <- total_no_computer_sj*100/total_sj_comp

This is the total number of households in San Jose without a computer:

print(total_no_computer_sj)
## [1] 18010

This is the percent of households in San Jose a computer:

print(perc_no_computer_sj)
## [1] 5.840027

This map highlights San Jose block groups by count of households without a computer.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_blockgroups_computer %>% 
      pull(`No computer`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_blockgroups_computer %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`No computer`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    label = paste0(sj_blockgroups_computer$`No computer`)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_blockgroups_computer$`No computer`,
    title = 'Households<br>Without a Computer'
  )

This is a comparable map, except that households without a computer are represented as a percent of the total number of households within the block group.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_blockgroups_computer %>% 
      pull(`percent_no_computer`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_blockgroups_computer %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`percent_no_computer`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    label = paste0(sj_blockgroups_computer$`percent_no_computer`)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_blockgroups_computer$`percent_no_computer`,
    title = 'Percent Without<br>a Computer'
  )

Households Below 80% AMI

Using area median income (AMI) is a more appropriate way to measure low-income status, as the cost of living in San Jose is substantially higher than the national average.

sj_median_income_by_block <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "B19013_001E"
  ) %>%
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  rename(
    Median_Income = B19013_001E 
  ) %>% 
  filter(!is.na(Median_Income)) %>%
  left_join(sj_blockgroups, by = c("blockgroup" = "GEOID")) %>%
  na.omit()
sj_ami_by_block <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*", 
    regionin = "state:06+county:085",
    vars = "group(B19001)"
  ) %>%
  mutate(
    blockgroup =
      paste0(state,county,tract,block_group)
  ) %>% 
  select_if(!names(.) %in% c("GEO_ID","state","county","tract","block_group","NAME")) %>% 
  dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
  group_by(blockgroup) %>% 
  summarize(
    Total = B19001_001E,
    `Under 75,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E),
    #sum(lapply(2:12, function(x) as.name(paste0("B19001_00",x,"E"))))
    `Under 100,000` = sum(B19001_002E, B19001_003E, B19001_004E, B19001_005E, B19001_006E, B19001_007E, B19001_008E, B19001_009E, B19001_010E, B19001_011E, B19001_012E, B19001_013E)
  ) %>% 
  mutate(
    `% under 75,000` = `Under 75,000` / Total * 100,
    `% under 100,000` = `Under 100,000` / Total * 100
  ) %>% 
  left_join(sj_median_income_by_block %>% dplyr::select(-Median_Income)
  ) %>%
  st_as_sf()

sj_ami_by_block[is.na(sj_ami_by_block)] = 0

total_under_ami = sum(sj_ami_by_block$`Under 100,000`)
total_sj_ami <- sum(sj_ami_by_block$Total)
perc_sj_under_ami <- total_under_ami*100/total_sj_ami

This is the total number of households in San Jose below 80% AMI:

print(total_under_ami)
## [1] 277191

This is the percent of households in San Jose without internet:

print(perc_sj_under_ami)
## [1] 43.61607

This map highlights San Jose block groups by count of households without a computer.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_ami_by_block %>% 
      pull(`Under 100,000`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_ami_by_block %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`Under 100,000`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    label = paste0(sj_ami_by_block$`Under 100,000`)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_ami_by_block$`Under 100,000`,
    title = 'Households<br>Under 80% AMI'
  )

This is a comparable map, except that households without a computer are represented as a percent of the total number of households within the block group.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_ami_by_block %>% 
      pull(`% under 100,000`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_ami_by_block %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`% under 100,000`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    label = paste0(sj_ami_by_block$`% under 100,000`)
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_ami_by_block$`% under 100,000`,
    title = 'Percent Under<br>80% AMI'
  )

K-12 Enrollment

This segment takes a look at the segment of the San Jose population 3 years and older enrolled in kindergarten through 12th grade.

sj_blockgroups_k12 <-
  getCensus(
    name = "acs/acs5",
    vintage = 2018,
    region = "block group:*",
    regionin = "state:06+county:085",
    vars = c("B14007_001E", "B14007_004E", "B14007_005E", "B14007_006E", "B14007_007E", "B14007_008E", "B14007_009E", "B14007_010E", "B14007_011E", "B14007_012E", "B14007_013E", "B14007_014E", "B14007_015E", "B14007_016E")
  ) %>%
  mutate(
    `School Enrollment` = B14007_004E + B14007_005E + B14007_006E + B14007_007E + B14007_008E + B14007_009E + B14007_010E + B14007_011E + B14007_012E + B14007_013E + B14007_014E + B14007_015E + B14007_016E
  ) %>%
  mutate(
    `Percent Enrollment` = `School Enrollment`*100/B14007_001E
  ) %>%
  rename(
    Total = B14007_001E,
    `Kindergarden` = B14007_004E,
    `Grade 1` = B14007_005E,
    `Grade 2` = B14007_006E,
    `Grade 3` = B14007_007E,
    `Grade 4` = B14007_008E,
    `Grade 5` = B14007_009E,
    `Grade 6` = B14007_010E,
    `Grade 7` = B14007_011E,
    `Grade 8` = B14007_012E,
    `Grade 9` = B14007_013E,
    `Grade 10` = B14007_014E,
    `Grade 11` = B14007_015E,
    `Grade 12` = B14007_016E,
  ) %>%
  mutate(
    blockgroup = paste0(state,county,tract,block_group)
  ) %>%
  #select_if(!names(.) %in% c("GEO_ID", "state", "county", "tract", "block_group", "NAME")) %>%
  #dplyr::select(-c(contains("EA"),contains("MA"),contains("M"))) %>%
  #gather(key = "variable", value = "estimate", -blockgroup) %>%
  #mutate(label = acs_vars$label[match(variable,acs_vars$name)]) %>%
  #dplyr::select(-variable) %>%
  left_join(sj_blockgroups, by = c ("blockgroup" = "GEOID")) %>%
  na.omit() %>%
  st_as_sf

total_k12 = sum(sj_blockgroups_k12$`School Enrollment`)
total_sj <- sum(sj_blockgroups_k12$`Total`)
perc_k12 <- total_k12*100/total_sj

This is the total number of children in San Jose enrolled in K-12:

print(total_k12)
## [1] 161712

This is the percent of San Jose population enrolled in K-12:

print(perc_k12)
## [1] 16.98577

This map highlights San Jose block groups by raw count of children in K-12. It is important to note that these numbers reflect 2018 estimates.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_blockgroups_k12 %>% 
      pull(`School Enrollment`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_blockgroups_k12 %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`School Enrollment`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    group = "K-12 Enrollment",
    label = paste0(sj_blockgroups_k12$`School Enrollment`)
  )  %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_blockgroups_k12$`School Enrollment`,
    title = 'Individuals Enrolled<br>in K-12'
  )

This is a comparable map, except that K-12 enrollment is represented as a percent of the total block group population.

pal <- 
  colorNumeric(
    palette = "Blues",
    domain = 
      c(0,sj_blockgroups_k12 %>% 
      pull(`Percent Enrollment`) %>% 
      unique())
  )

leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(
    data = sj_blockgroups_k12 %>%
    st_transform(4326),
    fill = TRUE,
    fillColor = ~pal(`Percent Enrollment`),
    color = "white",
    weight = 0.5,
    opacity = 0.5,
    fillOpacity = 0.5,
    group = "K-12 Enrollment",
    label = paste0(sj_blockgroups_k12$`Percent Enrollment`)
  )  %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = sj_blockgroups_k12$`Percent Enrollment`,
    title = 'Percent Enrolled<br>in K-12'
  )

Lack of Small Cell Coverage

sj_parcel_distance <-
  readRDS("sj_parcel_distance.rds")

# we want to start by assigning each parcel to a block group


# then we want to add up the average distance for each block group